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The low temperature thermal conductivity of a d-wave superconductor is considered for arbi- 
trary strength of the impurity potential. A random distribution of vortices is included within the 
usual semiclassical approach. The vortex-quasiparticle scattering is accounted for phenomenologi- 
cally, in addition to impurity scattering. Application of the theory to recent data on high purity 
YBa2Cu3 06.99 shows excellent agreement with experiment, both in zero field and in the presence of 
an applied magnetic field. We find that the strength of the impurity potential required to describe 
the data deviates qualitatively from either the unitary or Born limits. 
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Recent measurements of the in-plane thermal conduc- 
tivity (k) of a high purity YBa2Cu3 06.99 single crystal 
Q have revealed disagreement with the standard d-wave 
theory in several respects. The very low temperature (T) 
variation of k fits neither the unitary nor the Born hmit 
of impurity scattering theory. In addition, when a mag- 
netic field {H) is applied up to 13 Tesla, the T — data 
disagree with the predictions of the usual semiclassical 
theory of the Doppler shift which accounts for the super- 
currents circulating around the vortex cores, k initially 
increases very rapidly reaching almost twice its universal 
value kqo as H increases to a fraction of a Tesla. Beyond 
this value, k remains almost unchanged up to 13 Tesla 
while the semiclassical theory of Kiibert and Hirschfeld 
predicted it should continue to rise as y/H. 

Quasiparticles can be created not only thermally but 
also by the magnetic field. This additional creation of 
quasiparticles gives rise to a rapid increase of n. However, 
as the density of vortices is increased the quasiparticles 
should also experience additional vortex-quasiparticle 
scattering due to Andreev reflection, and this should 
have the corresponding effect of decreasing k. It is well- 
known that vortex scattering has important effects on 
the thermal conductivity of a type II superconductor 
0. However, it was not included in the work of Ref. 2], 
which applies only to the dirty limit, where it is assumed 
that impurity scattering is dominant. To treat the pure 
limit where the vortex scattering plays a significant role, 
we follow a semi-phenomenological approach. We use 
Matthiessen's rule to add the vortex scattering to the 
usual impurity term which is left unchanged. 

A random distribution of vortices is assumed and a 
configuration average is done as usual. After the averag- 
ing the vortex scattering is dependent only on the single 
vortex cross section. This quantity has been evaluated 
by Yu et al. 0], who gave a simple formula which is 
proportional to the magnetic energy Eh- Here the pro- 
portionality constant will be treated as a phenomenolog- 
ical parameter. Franz has also given a dimensional 
argument for this form. However, his approach does 



not give the correct behavior of k at low temperatures 
because in his calculation k — > kqq as T — > regard- 
less of H, which is inconsistent with experiments. While 
Vckhtcr and Houghton did include Andreev scatter- 
ing in their work, they start, however, with an Abrikosov 
lattice which introduces extra complications associated 
with Bloch's theorem. Their model does not apply to a 
random distribution and it also predicts a Vi? behavior 
at low temperature. Moreover, their theory is not well- 
justified for H ^ Hc2- Arguing that the Dirac nodes 
remain unaffected even in the vortex state (except for 
renormalized Fermi and gap velocities), Franz and Vafek 
concluded that for H^i ^ <C ffc2 the thermal con- 
ductivity is independent of H. This may happen only in 
a perfect vortex lattice; however, it is debatable if vor- 
tices form a lattice in a rigorous sense even in a very pure 
sample. 

In the work of Hill et al. 1] for H ~ 0, k/T rises 
out of its universal limit very fast as a function of tem- 
perature, and reaches 4 times this universal value within 
AOOmK . This behavior cannot be understood within ei- 
ther of the two limits extensively studied in the past 
Born and unitary limits, corresponding to the impurity 
potential [/; — > and J7i — > oo, respectively 0. In this 
paper we consider explicitly the intermediate case. A 
dimensionless parameter c is introduced, as usual, with 

— l/{UiN{Oj), where N{0) is the normal state density 
of states at the Fermi level. Within a i-matrix approx- 
imation, the effective scattering in the superconducting 
state can be characterized by c and F — ni/(7riV(0)), 
where rii is the impurity density. We find that the data 
of Hill et al. 3| can be fit naturally with a value of 
c ~ 0.014, which does not belong to either of the two 
limiting cases. The rapid rise in k/T with T also reflects 
the high purity of the sample {i.e small value of F). 

An important feature of the impurity scattering is that 
the scattering rate 7(0;) in the superconducting state is 
strongly dependent on frequency and that our fit to data 
for YBaaCuaOe.gg with ~ 90K gives 7(0) = 0.29K, 
which is very small in comparison to, for example, the 
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value of the magnetic energy at Hd- This value is of 
order 10"^Ao while 7(0) ~ O.OOlSAo, where Aq is the 
amplitude of a d-wave gap and is 2.14Tc for a d-wave 
superconductor. Thus when we include a magnetic field 
H into the theory we will need to work in the regime 
Eh > 7(0)- The limit Eh < 7(0) is unreahstic for such 
pure samples. Nevertheless, provided Eh/'^o and T /An 
remain small, we can apply the nodal approximation [iflj 
to reduce complications associated with the calculations. 

The superconducting (2 x 2) matrix Green's func- 
tion in the Nambu notation takes the form G(k, = 
(wfo + AkTi + ^k'r3)/(w^ - Ck ~ ^k)' where f's are the 
Pauli matrices in spin space, Ak a d-wave order param- 
eter, and ^k is the electronic energy dispersion in the 
normal state. The renormalized Matsubara frequencies 
Lu{iLUn) ~ iiUn — S(d)) includc the self-energy S(ti)), which 
depends on what interactions are taken into account. 
Here we consider only the scalar component of the self- 
energy assuming that other components can be neglected 
or absorbed into ^k or Ak^3|- our consideration for 
a vortex state, it has two parts: the impurity part 
and the vortex part T.y . For we use the self-consistent 
i-matrix approximation while for E^, we take a more phe- 
nomenological approach due to the lack of a detailed mi- 
croscopic footing. 

In the absence of a magnetic field, lu(u! + i6) — uj — 
Yii^ret{'^)- Introducing the spectral functions for the 
diagonal [G(k, iw)] and off-diagonal [i^(k, zlj)] compo- 
nents of G(k, iw), A(k,u;) = -2Im [G(k, w -I- and 
i3(k, o)) = — 2Im [i^(k, a; -I- respectively, we obtain 
the thermal conductivity at T 
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where ^/(g) is the Fermi (gap) velocity and /(w) is the 
Fermi function. As we mentioned earlier, application of 
the nodal approximation allows us to carry out the sum- 
mation over k and yields 
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arctan[/i(cj)]. Here h{uj) 



where A{tj) = 1 + h{uj) + 

is defined as h{uj) = {uj — Re [Ei_ret])/7(^) with 7(w) = 
— Im [E,; ret]- In order to calculate the thermal conductiv- 
ity in Eq. |2Jl, we need to determine the self-energy Ei^^et 
self-consistently. In the i-matrix approximation E^^^et 
can be written as Ei_ret('i) = rGo((I')/ [c^ — Go(tl')] , 
where Go = [2TrN{0)]-^ Xlk Tr[foG^et]- As we described, 
c here controls the strength of the impurity scattering. 
Applying the nodal approximation, we obtain 
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Substituting Eq. ^ into Ei^^et, we compute Re[Ei^ret] 
and Im[Ei_j.et] self-consistently. In the presence of a mag- 
netic field, the results for Ej ret remain unchanged. 

Since A{uj) 2 as T — > 0, it is easy to 

show that k{T)/T approaches its universal value, 
(2/3) [vf /vg + Vg/vf] = Koo/T. In Fig. 1 we show re- 
sults for k{T)/T normalized by the universal value for 
different values of c. The figure is motivated by the data 
of Hill et al. on a very pure sample of YBa2Cu3 06.99. 
Their experimental fit appears as solid circles correspond- 
ing to (0.16 + 3.0r2)/0.16 Q. We should note that this 
data requires subtraction of a phonon contribution; this 
can be done in a variety of sensible ways, which leads 
to an uncertainty in the electronic part of several per- 
cent The solid curve shown, a theoretical calcula- 
tion with c = 0.014, is thus in excellent agreement with 
the data. While there are two parameters in the impu- 
rity theory, F and c, a constraint was applied to all fits: 
namely, 7(0) is set to be Q.29K . This constraint comes 
from consideration of finite as well as zero magnetic field 
data as will be explained later. As seen in Fig. 1 the uni- 
tary limit does not fit the data at finite T. On the other 
hand we can achieve reasonable agreement with the data 
for c = 0.014 and F/Aq ~ 5.05 x 10^''^ over the entire 
temperature range shown in Fig. 1. This corresponds to 
a very clean sample. If F were to be interpreted as the 
normal state impurity scattering rate, the corresponding 
mean free path would be anomalously large. We take this 
as an indication that for this sample the underlying nor- 
mal state, on which the superconductivity is built, may 
not be an ordinary Fermi liquid. 

We note that the theory does not give an exact straight 
line while the experimental data has been parameterized 
with a linear fit P|. However, it should be remembered 
that to obtain the electronic thermal conductivity from 
experiment it was necessary to subtract a considerable 
phonon contribution. Uncertainty remains about this 
subtraction as we will elaborate on later and so the small 
disagreement that remains between theory and experi- 
ment is not a serious limitation on the fit obtained. In 
the inset of Fig. 1 we show 7(0;) obtained for c = 0.014 
(solid curve) compared with the unitary limit (dotted 
curve). The curves are completely different. They both 
start with 7(0)/Ao = 0.0015 by design but for finite lo 
the solid curve rises while the dotted curve drops mono- 
tonically. These differences lead to the large difference in 
the T dependence of k{T)/T shown in Fig. 1. 

In the presence of a magnetic field. Next we wish to 
include the effect on the thermal conductivity of an ex- 
ternal magnetic field H applied along the c-axis. To this 
end we need to know the self-energy of a quasiparticle 
due to vortices E^, as well as for the impurities. For- 
mally speaking, the semiclassical approximation takes 
into account an approximate real part of the self-energy: 
Re [E„] ~ Vs(r) • k where Vs(r) is the superfluid veloc- 
ity at a position r measured from the center of a vortex 




FIG. 1: The thermal conductivity n{T)/T in units of its uni- 
versal value as a function of (T/Tc)^ in the very low T <^Tc 
range. In aU cases, Tc = and 7(0) = Q.2'dK. The impu- 
rity concentration parameter V varies as c is changed. Each 
curve is labeled by the c value used. The experimental fitting 
appears as solid circles. The inset shows the self-consistent 
results for ^{u)) for c = (dotted) and c = 0.014 (solid curve). 



core. Here r is restricted to the vortex unit cell and 
an average over r is to be carried out. For the imagi- 
nary part of Yiy, we consider the Andreev reflection as a 
possible mechanism. When a quasiparticle moves from 
r to r', the Andreev reflection takes place effectively if 
+ + v,(r) • k = |Ak| + Vs(r') • k. It has been 
shown Q that for nodal quasiparticles the vortex scat- 
tering rate 71, is determined mainly by the intervortex 
distance R; namely, 7„ ~ Vf/R. This implies that within 
the nodal approximation a vortex plays the role of a de- 
fect or a scattering center. Since the intervortex distance 
R ~ y/^a/H, where $0 is a flux quantum, and the mag- 
netic energy Eh is given by Eh = aJ t: /2v f H /'^i^^ 
where a is a parameter of order unity 2] , we can write 
7„ = hEn^ which defines a vortex scattering parameter 
h. Instead of evaluating the microscopic quantities which 
determine the size of b from first principles, we propose 
to treat it as a phenomenological parameter which we 
can use to fit the zero temperature data of the thermal 
conductivity as a function of H . 

An averaged value of the thermal conductivity over the 
vortex distribution can be evaluated as usual: n{H^ T) ~ 
J deV{e) K{e,T), where V{e) = / dr 5 [e - Vs(r) • k] 
is the vortex distribution within the area of a vortex unit 
cell A. The appropriate average (series or parallel) de- 
pends on the direction of the field as discussed by Kiibert 
and HirschfeldQ. We will not deal directly with these 
issues here. We made numerical calculations assuming 
series and parallel arrangements for the thermal conduc- 
tivity but have found no qualitative differences. Now the 



FIG. 2: The normalized k{T)/T as a function of H for the 
BZO sample (triangle) and the YSZ sample (circle) . The inset 
gives the variation of the magnetic field saturation function 
gib) vs. b. (See text for more explanation). 



thermal conductivity in a vortex state becomes 
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where A{lo, e) = 1 -I- /i(ti;, e) -I- 



arctan[/i(Li;, e)]. 
The effects 



Note the similarity between Eqs. Q and Q 
of the Doppler shift and the Andreev scattering appears 
in h{uj,e); h{uj,e) = (cj - Re [Si^^et] - e) /7t(^), where 
7t(a;) = j{uj) + bEn. 

Taking the limit of T = in Eq. (@J leads to a simplified 
expression as follows: k{H,0) = (koo/2) / deV{€)A{0,e) 
with ^(0, e) = 1-f [e/7t(0) + 7t(0)/e] arctan [e/7t(0)]. To 
proceed further some vortex distribution function 'P{e) 
needs to be specified. We consider the Gaussian distribu- 
tion; V{e) = exp [-e'^/Ej^] /{^Eh), which is believed 
to be favorable for Andreev scattering Now we ob- 
tain the thermal conductivity as T — > 0: 
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where h{0,x) = x/ [yiO)/ Eh + b]. For Eh > 7(0), 
h(0,x) ~ x/b and, thus, k(0,H)/koo ~ 5(^)1 where 
g{b) = (l/2v^) / dxe"^' [1 + {x/b + b/x) arctan(a;/6)]. 
Consequently, the magnetic field saturated value of 
k{0,H) is g{b)Koo and not kqo as found in Ref.j^. 

The function g{b) is shown in the inset of Fig. 2. Note 
that as 6 ^ 0, g{b) goes to infinity as y/7r/{4b) and no 



saturation occurs in this case. This also corresponds to 
the dirty hmit case considered in Ref.0| where it is as- 
sumed that 7(0) is much bigger than bEn- In the main 
frame of Fig. 2 we show our results for the normalized 
k{H, T)/T vs 7? at T = for the two samples considered 
by Hill et al. , a very pure sample denoted by BZO with 
Tc ~ 90K (triangle) and an impure sample denoted by 
YSZ (circle). In both cases the agreement with the data 
is excellent (compare with Fig. 3 of Ref. 0]). 

For the BZO sample with 7(0) ~ 0.29K the value of b 
obtained in the fit is 6 = 0.3 assuming w/ — lO^cra/s and 
a 1/2 so that Eh /Aq ~ 0.036/ff T-^/^Jll. Using 
Hci ~ 30mT of YBa2Cu30:r with Tc ~ 91K [ij we have 
Eh/Ao ~ O.OO6A0, which is greater than 7(0). For the 
dirty sample (YSZ) with Tc ~ 93K, the estimated value 
of 7(0) is about 13 K ^. Using these parameters with 
~ 2 X lO^cm/s, we found that no vortex scattering 
parameter needs to be introduced to fit the data; namely, 
6 = 0. In this case the magnetic energy at a few Tesla is 
still of the order of 7(6) and so the data is not sensitive 
to a small value of b. 

In Fig. 3 we show results for the normalized k{H, T)/T 
in units of the universal value as a function of {T/TcY ■ 
The solid curve reproduces our H — result for c — 
0.014. The different symbols represent different values of 
H; namely, if = IT (open circle), 4T (open square), and 
9T (open diamond). Note that within 10% or less the 
saturated value of k{H,T)/T due to H is temperature 
independent. This justifies the process used in Hill et 
al. to subtract out the phonon contribution from their 
raw data which includes both an electronic and a phonon 
part. We also emphasize that the crossing with the H = 
curve occurs around {T/Tc)^ = 5 x 10"^ or T ~ 0.2/^, 
which agrees well with experiment. We point out that to 
get this agreement we insisted on 7(0) ~ 0.29i^ and the 
crossing temperature T* ~ 0.2K is close to 7(0) so that 
it is the crossing that determines this parameter in our 
approach. For 6 = the crossing cannot occur. This im- 
plies that both the crossing itself and the saturated value 
of k(0, T) show the significance of the vortex scattering. 
In the inset we show a few more results for small values 
of Hci). This inset shows that for a weak field the 
thermal conductivity retains significant temperature de- 
pendence although considerably reduced from the H = 
case. As one can see in this inset or the main frame 
of Fig. 2, at T = the thermal conductivity depends 
strongly on H for Hd ^ H < IT. This can be seen in 
the experiments. 

In summary we have given a theory of the low temper- 
ature thermal conductivity of a d-wave superconductor 
using a self-consistent treatment of impurity scattering 
intermediate between the unitary and Born limits. We 
have also included the effects of a magnetic field in the 
range Hd < H ^ Hc2 within the semiclassical approx- 
imation which takes into account the Doppler shift due 
to the circulating supercurrent around the vortex cores. 
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FIG. 3: The normalized n{T)/T as a function of (T/Tcf. 
The solid curve is our result for H = with c = 0.014 repro- 
duced from Fig. 1. The open circle, square, and diamond are 
for H — 1,4, and 9T, respectively. The inset shows additional 
results for H = 0.04T, 0.16T as well as ii" = and = IT 
for comparison. 



Assuming a random distribution of vortices, the Andreev 
scattering from the vortices is considered and this leads 
to saturation in the thermal conductivity at T = for 
a sufficiently large value of H. We find good qualitative 
and quantitative agreement in all respects with recent 
data on a very pure sample of YBa2Cu3 06.99. We also 
verify that the phonon contribution can be subtracted 
out using the thermal conductivity for a sufficiently high 
magnetic field. 
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